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ABSTRACT 


The  large  amplitude  vibrations  of  thin  elastic  plates  and  shallow 
shells  having  various  boundary  conditions  and  subjected  to  random  exci¬ 
tation  are  investigated  by  using  various  approximate  techniques. 

The  random  vibrations  of  rectangular  plates  and  circular  plates 
subjected  to  white  random  excitation  are  simulated  numerically  by  two 
different  methods.  The  first  method  is  that  the  governing  equations 
are  reduced  to  a  single-degree-of- freedom  dynamical  system  and  the 
reduced  equation  is  then  integrated  numerically  by  the  Runge-Kutta 
method  employing  the  simulated  approximate  white  noise  as  an  input. 

The  second  method  consists  in  integrating  the  equation  of  motion  and 
the  compatibility  equation  numerically  by  a  finite-difference  method 
employing  the  simulated  approximate  white  noise  as  an  input.  To  compare 
the  results  obtained  by  the  simulation  methods  with  those  by  other  methods, 
the  single-degree-of- freedom  system  equation  is  solved  exactly  using  the 
Fokker-Planck  equation,  and  solved  approximately  by  the  equivalent  lineari¬ 
zation  technique.  Also  presented  is  the  response  analyses  of  shallow  shell 
to  white  noise  by  (1)  numerical  simulation  using  the  single-degree-of- 
freedom  equation  and  (2)  the  Fokker-Pl anck  equation.  It  has  been  shown 
that  the  solutions  by  the  numerical  simulation  are  close  to  those  obtained 
by  the  equivalent  linearization  technique  and  the  Fokker-Planck  approach 
while  the  second  numerical  simulation  gives  rather  poor  solutions. 
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INTRODUCTION 


In  many  physical  problem  areas  there  are  situations  where  a  mechanical 
system  is  excited  by  a  random  type  load  and  the  response  of  the  system 
displays  a  random  trend.  Since  the  response  process  is  not  deterministic, 
the  response  analysis  must  be  treated  statistically. 

During  the  last  two  decades  much  research  effort  in  the  area  of  vibrations 
has  been  devoted  to  the  investigation  of  the  structural  response  to  random 
excitation.  Motivation  for  such  research  has  arisen  due  to  the  development 
of  large  jet  engines  and  rocket  motors  which  produce  random  pressure  fields 
of  high  intensity.  Since  the  level  of  random  excitation  generated  by  jet 
aircraft  and  missiles  provides  a  severe  environment  with  respect  to  fatigue 
failure  of  structures,  the  investigation  of  the  response  of  structures  to 
random  excitation  plays  an  important  role  in  the  fields  of  aircraft  and 
missile  design. 

The  works  of  Crandall  [1,2],  Bolotin  [3],  Crandall  and  Mark  [4],  and 
Lin  [5],  contain  various  topics  in  random  vibration  analysis.  So  far  the 
literature  in  the  area  of  random  vibration  analysis.  So  far  the  literature 
in  the  area  of  random  vibration  has  been  surveyed  by  Crandall  [6],  Smith 
[7],  Bolotin  [8]  and  Vorovich  [9].  Therefore  in  the  present  Chapter  only 
random  vibration  of  thin  elastic  plates  and  shallow  shells  will  be  reviewed 
briefly. 

Analysis  of  the  response  of  a  linear  elastic  structure  with  known  normal 
modes  to  random  excitation  of  a  given  power  spectrum  have  been  carried  out 
by  many  investigators  [10  through  13].  Usually  the  input  random  load  is 
assumed  to  be  stationary,  ergodic,  Gaussian  with  zero  mean  value.  For  such 
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a  case  for  a  linear  system  it  is  possible  to  relate  the  statistical  des¬ 
cription  of  the  output  to  that  of  the  input.  Hence,  the  mean-square 
response  can  be  evaluated  under  the  assumption  of  a  stationary  Gaussian 
input.  However,  the  results  of  such  analysis  are  valid  only  for  small 
lateral  displacements.  For  moderately  large  displacements  it  is  neces¬ 
sary  to  take  the  effects  of  nonlinearities  into  account. 

The  motion  of  geometrically  nonlinear  elastic  thin  plates  and  shallow 
shells  is  described  by  a  system  of  two  coupled  nonlinear  partial  differen¬ 
tial  equations  in  terms  of  the  lateral  displacement  w  and  the  stress 
function  F,  the  so-called  dynamic  analog  of  the  von  Karman  equations.  The 
load  is  again  considered  to  be  a  stationary  Gaussian  random  process  with 
zero  mean  value.  Unfortunately,  no  exact  solution  for  this  problem  has 
been  found.  Only  approximate  solutions  are  possible.  One  approximate 
solution  is  to  reduce  the  partial  differential  equations  to  a  system  of 
ordinary  differential  equations  for  the  generalized  coordinates  and  to 
obtain  an  approximate  solution  by  employing  techniques  used  in  nonlinear 
mechanics.  To  do  this,  assuming  that  the  lateral  displacement  w  is 
expressed  as 


N 

w=  Z  f,(t)T.(x,y) 
i  =  l  1  1 


(1-1) 


where  f. (t)  are  generalized  coordinates  and  T^(x,y)  is  the  coordinate 
function  representing  the  deflected  shape  of  the  structure,  and  using 
approximation  technique  (for  example  Galerkin's  method)  the  problem  is 
reduced  to  a  system  of  N  dynamical  equations  with  respect  to  f^(t),  i.e.. 


(1-2) 


N 
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where  c  =  generalized  damping  term 

g.j  =  generalized  nonlinear  stiffness  term 
Qi  =  generalized  random  forcing  function 
An  alternative  form  of  the  governing  equations  for  f.  (t)  derived  by 
mean  of  the  Lagrangian  equation  takes  the  form 


d'fi<t)  .  +_9U_ 

dt2  dt  3-f. 


Qi(t) 

=  1,2,  ...,N 


(1-3) 


where  U  is  the  potential  energy  of  the  system. 

Only  under  the  restriction  that  the  generalized  random  forcing  function 
Qn*  (t )  is  stationary,  Gaussian,  white  noise  or  filtered  white  noise,  can 
eqs.  (1-2)  or  (1-3)  be  solved  exactly  in  terms  of  the  joint  probability 
density  using  the  Markov  process  and  the  associated  Fokker-Pl anck  equa¬ 
tion  [14].  This  approach  was  first  applied  to  the  case  of  a  nonlinear 
system  by  Andronov  et  al  [15].  Herbert  [16]  investigated  the  multi - 
mode  response  of  beams  and  plates  to  white  noise  excitation  using  this 
approach  and  showed  that  the  probability  density  function  of  the  model 
amplitude  is  non-Gaussian  and  statistically  dependent.  Dimentberg  [23] 
applied  the  approach  to  the  curved  panel  problem  taking  one  term  of 
the  series  for  normal  deflection  in  eq.  (1-1)  and  studied  the  fatigue 
damage. 

If  the  random  forcing  is  not  assumed  to  be  white  noise,  only  an 
approximate  solution  to  eq.(l-2)  or  (1-3)  is  possible.  One  approximate 
technique  for  this  type  of  equation  is  the  equivalent  linearization 
technique  which  was  originated  by  Krylov-Bogol iubov  [17]  in  determinis¬ 
tic  theory  and  was  applied  to  problems  of  random  vibrations  by  Booton 
[18]  and  Caughey  [19].  Lin  [20]  investigated  the  single-mode  response 
of  a  flat,  plate  undergoing  moderately  large  deflections  subject  to 
stationary  Gaussian  excitation  whose  power  spectrum  is  relatively  flat. 
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Another  technique  to  obtain  an  approximate  solution  to  eq.(l-2) 
is  the  perturbation  method.  This  technique  was  first  introduced  to 
random  vibration  problems  by  Crandall  [21].  This  approach  can  be 
applied  only  to  systems  with  very  small  nonlinearities. 

Besides  the  above  mentioned  approximate  approach,  another  method 
to  estimate  the  response  of  a  nonlinear  elastic  structure  to  random 
excitation  is  a  numerical  simulation  technique  which  has  been  exten¬ 
sively  used  in  the  investigation  of  structural  response  due  to  earth¬ 
quakes.  This  technique  is  to  digitally  simulate  a  physically  realiza¬ 
ble  random  load  and  to  integrate  the  equation  of  motion  numerically 
by  employing  the  simulated  random  load  as  a  forcing  function.  Then 
we  compute  the  desired  statistical  properties  of  the  response.  Belz 
[22]  used  this  technique  to  investigate  the  problem  of  a  beam  subjected 
to  a  concentrated  random  driving  force.  However,  this  method  is  time 
consuming  if  the  structural  model  is  complex.  Instead  of  solving  the 
original  partial  differential  equation,  eq.(l-2)  may  be  integrated 
numerical  ly. 

In  the  present  study,  the  large  amplitude  vibrations  of  plates 
and  shallow  shells  having  various  boundary  conditions  and  subjected 
to  white  random  excitation  will  be  investigated  by  using  various 
approximate  methods.  Before  going  on  to  the  analyses  of  the  response 
to  random  excitation,  an  investigation  of  the  structural  response  to  a 
deterministic  load  will  be  made  as  a  preliminary  study  in  Chapter  II. 

In  Chapter  III,  a  Gaussian  stationary  white  random  process  is 
digitally  simulated. 

In  the  second  section  of  Chapter  IV,  the  random  vibrations  of 
rectangular  plates  and  circular  plates  subject  to  white  noise  are 
simulated  numerically.  Two  different  simulations  are  presented: 
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(1 )  the  governing  equations  are  reduced  to  a  single-degree-of- 
freedom  dynamical  system  by  taking  one  term  of  the  series  representing 
the  normal  deflection  in  eq.(l-l)  and  the  equation  corresponding  to 
the  case  of  N-l  in  eq.  (1-2)  is  then  integrated  numerically  by  the 
Runge-Kutta  method  employing  the  simulated  white  noise  as  an  input. 

(2)  the  equation  is  integrated  numerically  by  a  finite-difference 
method  employing  the  simulated  white  noise  as  an  input.  To  compare 
the  results  obtained  by  the  simulation  method  with  those  found  by 
other  methods,  the  single-degree-of  freedom  system  equation  is  solved 
exactly  using  the  Fokker-Planck  equation.  Also  the  approximate  solu¬ 
tions  obtained  by  the  equivalent  linearization  technique  are  presented. 
In  the  third  section  of  Chapter  IV,  the  response  analysis  of  shallow 
shells  to  white  noise  is  carried  out  by  (1)  numerical  simulation  using 
the  singl e-degree-of-freedom  system  equation  and  (2)  the  Fokker-Planck 
equation. 
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II 

RESPONSE  OF  THIN  ELASTIC  PLATES  AND 
SHALLOW  SHELLS  TO  STEP  LOAD 

In  this  Chapter,  the  response  analysis  of  flat  rectangular  plates, 
flat  circular  plates,  aribtrary  shallow  shells  with  rectangular  boun¬ 
daries  and  shallow  spherical  shells  with  a  circular  boundary  to  a  uni¬ 
formly  distributed  step  load  will  be  discussed.  Nonlinear  partial 
differential  equations  governing  the  finite  amplitude  deflections  of 
plates  and  shallow  shells  are  approximated  by  the  finite-difference 
equations  by  use  of  the  Crank-Nicol son  finite-difference  scheme  [24] 
and  these  difference  equations  are  then  solved  numerically  using  a  CDC 
3600  digital  computer.  The  computer  program  for  this  analysis  can  be  used 
for  investigation  of  response  to  an  arbitrary  input  forcing  function  if 
the  input  forcing  function  is  digitally  simulated.  In  Chapter  IV  the 
computer  program  written  in  this  study  will  be  used  for  digital  simulation 
of  random  vibrations  of  plates.  The  purpose  of  the  study  in  the  present 
Chapter  is  that  before  simulating  random  vibrations  we  determine  the 
numerical  stability  of  the  solution  and  the  accuracy  of  the  solution  by 
investigating  the  response  to  a  step  load.  Since  the  steady  state  response 
of  the  damped  system  to  step  load  must  agree  with  that  under  static  load, 
we  can  check  the  computer  program  and  the  accuracy  of  the  approximate  solu¬ 
tions  by  comparing  the  values  obtained  here  with  existing  results  for  the 
static  load. 

2.1  Rectangular  Plates  and  Shallow  Shells  of  Rectangular  Contour 

2.1.1  Analysis 


Governing  Equations 

Consider  a  thin,  elastic,  isotropic  shallow  shell  rectangular  in  plan 
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with  double  curvatures  and  of  constant  thickness  h.  Young's  modulus  E, 
and  Poisson's  ratio  v.  The  origin  o  of  the  curvilinear  coordinates  x-y-z 
is  chosen  at  a  point  of  the  middle  surface  corresponding  to  one  of  the 
corners  of  the  shell.  (See  Figure  2-1)  Let  the  oz  axis  extend  along 
the  normal  to  the  middle  surface  toward  the  center  of  curvature.  The  ox 
and  oy  axes  are  drawn  parallel  to  the  lines  of  principal  curvature  of  the 
shell,  a  and  b  denote  the  dimensions  of  the  shell  along  the  ox  and  oy 

axes.  Also  k  and  k  are  the  curvatures  of  the  shell  which  remain  constant 

x  y  * 
along  the  ox  and  the  oy  axes,  respectively.  Let  w  be  the  displacement 

of  a  point  in  the  middle  surface  along  the  oz  axis. 

The  differential  equations  governing  the  finite  amplitude  vibrations 

of  such  a  shell  are  [25] 


ph 


3*w  +  c*ph|£—  +  DV4w 


3F" 


h32F*  32w*  s 

h“3yT  (kx  lx7-) 


.  ha2F*/k  +  2h 92f*  aV  * 

+  h1)F^ky  ly1^  "  2h3x3y  3x3y  q 


(2-1) 


rjbr*  _  rr/32W*N2  .  k 

v  h  tLN5x5y'  3x2  3y2  x  3y2 


.  ★  .  ★ 

-  k..^r] 


y  3x2- 


(2-2) 


where  D=Eh 3/l 2 ( 1  -v2 )  is  the  flexural  rigidity  of  the  shell,  F  is  the  Airy 

★ 

stress  function,  p  is  mass  density  of  the  material,  c  is  the  damping  co- 

★ 

efficient  (assumed  to  be  constant),  q  is  the  lateral  load  and  t  denotes 
time. 

★  ★  ★ 

Membrane  stresses  a  ,  a  and  t  in  the  middle  surface  are  given  by 

x  y  xy 
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*  _  a2F 

°x  ~  3y2 

*  a2F* 

Gy  = 


3x2 


(2-3) 


xy 


32F 

3x3y 


The  strains  in  the  middle  surface  are  expressed  in  terms  of  F  as 


*  1  /32F 

ex  E  3y2 


32F  x 
"  V  3x2 


_  ★  _  ★ 

*  _  1  ,32F  32F  x 

ey  E'  3x2  "v  3y2 


(2-4) 


*  =  .  2 ( 1  +v )  12_F_ 


'xy 


3xay 


Considering  u  and  v  which  are  the  displacements  in  the  middle  surface 

in  the  x  and  y  directions,  respectively,  the  strains  are  expressed  in 
*  *  * 

terms  of  u  ,  v  and  w  as 


*  3u 
£  = 


3x 


kxw*  +^2 


*  i5L.  _  k  w*  +  !(^)2 
y  3y  y  2l3y  ] 


(2-5) 


*  3u*  3v*  3w*  3w* 

Yv  =  - -  +  - -  + 


"xy  3y  3x 


3x  3y 


The  bending  moments  are 

h  _  r>  /  3 2  w  ,  32w  \ 
i1x"'D(^P_+  Ty~) 


(2-6) 


My  -  -D( 


w 


9y2 


9  * 

9  w  x 
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Non-dimensional i zed  Equations 

Let  us  introduce  the  following  dimensionless  parameters: 


Or 


-  ^ 


K0  -  12(1  -v2 ) 


Txy  E  V 


n 


=  iL 
b 


u 


u  b 


T 


t(-T-)2 

phb4 


v 


v  b 


c 


C*(^ 


4  1 

)7 


w 


F  = 


EhT 


„  _  X/bv2 

°i  "  r-(h} 


*  4  4 

q  =  q  b^/Eh^ 

X  =  a/b  =  aspect  ratio 
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(2-6) 

cont'd 


(2-7) 


Using  the  above  dimensionless  parameters,  the  equation  of  motion  (2-1)  and 
the  compatibility  equation  (2-2)  are  now  expressed  by 
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a2w  .  aw 
at2  c  at 


4  4 

f3  w  ,  o_3 _ 

laf3T  ^2an2 


u  a2F/32w 


+ 


+  K 


a2F 

'oac 


/ifw  + 

2'an2n 


k2)  - 


2K 


a2F  a 


2w 


‘o a ^ 3 n  a^ari 


(2-8) 


a4F  +  2  a4F  +  a^F 

a?4  a?23n2  an4 


(3fw_  ^2  _  afw  a^_ 

a^an  as2  an2 


a2w 

■KiaTF 


.  a2w 

'2a|^ 


(2-9) 


Boundary  Conditions  and  Initial  Conditions 

In  this  study,  shallow  rectangular  shells  with  the  following  two 
boundary  conditions  are  considered: 

a)  All  four  edges  are  simply  supported  and  immovable  constrained 
against  in-plane  translation. 

b)  All  four  edges  are  clamped  and  immovably  constrained. 

Let  us  formulate  each  set  of  boundary  conditions.  For  the  simply  supported 
case,  the  deflection  w  along  four  edges  must  be  zero  and  there  is  no  bend¬ 
ing  moment  along  any  edge.  Thus,  the  analytical  formulations  of  these 

boundary  conditions  in  dimensionless  form  are  from  eq.(2-6), 

a2w 


w  =  0  and 


an2 

a2w 

a?2 


a2w 

X2’ 

a2w 

an2 


=  0 


=  0 


at  n  =  0  and  n  =  1 


w  =  0  and 


at  £  =  0  and  £  =  X 
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Since  w=0  along  n=0  and  n=l »  must  be  zer0-  Also, 

becomes  zero  along  5=0  and  5=X.  The  above  conditions  can  therefore 

be  written  as 

w  =  0  and  =0  at  n  =  0  and  n  =  1 

(2-10) 

w  =  0  and  =0  at  5  =  0  and  5  =  X 

If  the  edges  are  also  immovably  constrained  at  the  supports,  the  normal 
strain  in  the  middle  surface  parallel  to  the  edge  must  be  zero  along  the 
edge.  The  boundary  conditions  are  from  eq.(2-4). 


92F 

952 

92F 

V952" 

0 

at 

5=  o 

and 

92F 

9n2 

d2r 

w 

0 

at 

o 

II 

CT 

and 

zi 

n 

One  additional  condition  required  is  that  the  relative  displacement  of  the 
points  on  edges  5=0  and  5=X  for  any  given  value  n  is  equal  to  zero. 


u 


S=A 


rX 
■’  o 


(ff)  «  ”  0 

^  n=constant 


(2-12) 


From  eqs.(2-4)  and  (2-5),  we  have 


9 u  _  92F  „  3*F  l/9Ws2  ,  kl  w 
95  9ri2  ~V  af’2  "  ?  ar'  T? w 


2'95' 


Substituting  this  expression  into  eq . (2-12)  gives 
0  '  v  ^  -  - 

0 


f  0  -  fr)  *  ■  I  Ci(|f)2-^  W]  d5 


n  -  const. 


n  =  const. 


(2-13) 
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In  a  similar  fashion,  we  have 


i 

o 


/32F  d2Fs 


dn  = 
E,  =  const. 


[^■(|^)2  -  k2w]  dn 

5  =  const. 


Next,  consider  the  boundary  conditions  for  clamped  shells  with 
immovably  edges.  In  this  case  the  deflections  along  the  boundary  cr 
zero  and  the  plane  tangent  to  the  deflected  middle  surface  does  not 
at  the  edges. 

Therefore,  we  have 


w  =  0 


and  |*=o 

3n 


at  n  =  0  and  n  =  1 


w  =  0  and  =  0  at  £  =  0  and  £  =  A 

The  conditions  (2-11),  (2-13)  and  (2-14)  must  also  be  satisfied. 
For  initial  conditions  we  assume  the  body  to  be  at  rest. 


Then 


w  =  0 


and  |^=0 

3t 


at  x=  0 


The  transverse  load  applied  to  the  shallow  shell  in  this  Chapter 
a  uniformly  distributed  step  function  type  load  which  is  expressed  b; 
q(£,n,x)  =  Q  u  (t) 

where 


Q 

u(x) 

u(t) 


=  amplitude  of  the  load 
=  unit  step  function  defined 
1  when  x  >  0 


2  when  x  =  0 


by 


(2-14) 

e 

rotate 


(2-15) 


(2-16) 

is 

(2-17) 


(2-18) 


0  when  x  <  0 
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The  problem  consists  in  determining  the  functions  w  and  F  which 
satisfy  eqs.(2-8)  and  (2-9)  together  with  the  prescribed  boundary  con¬ 
ditions  and  initial  conditions.  Because  of  difficulty  in  solving  these 
simultaneous  partial  differential  equations  analytically,  these  equations 
will  be  solved  numerically  by  a  finite-difference  method.  The  idea  behind 
numerical  integration  of  eq.(2-8)  by  the  finite-difference  method  is  the 
following:  If  the  deflection  w  in  the  middle  surface  is  specified  at  a 
certain  time  level  x^,  the  stress  function  F  and  the  membrane  stresses 
are  determined  by  the  compatibility  equation  (2-9)  and  boundary  conditions. 
Given  membrane  stresses  at  time  level  x^,  we  seek  the  deflection  w  at  the 
next  level  T (<+ATwi th  load  Q  under  the  assumption  that  the  membrane  stresses 
at  each  point  remain  constant  in  the  time  interval  [t^^t^+At)].  By  this 
assumption,  the  equation  of  motion  (2-8)  is  treated  as  a  linear  partial 
differential  equation  in  the  time  interval  [x^,  (x^+Ax)]. 

Equation  (2-8)  is  reduced  to  two  lower  order  differential  equations 
by  introducing  the  two  variable  W  and  V  defined  as  follows: 


W  =  3fw 

H2  8n2 


(2-19) 


V  = 


3w 

3x 


(2-20) 


If  we  differentiate  eq.(2-19)  with  respect  to  x  and  substitute  eq. 
(2-20) ,  we  have 


3W  _  32V  3 2V 

HFr  ■572* 

which  is  the  first  desired  equation. 


(2-21) 

Substituting  eqs.(2-19)  and  (2-20) 


into  eq.  (2-9)  to  get. 


3V  +  c  V 
3x  cv 


-  ( 


32W 

3£2 


+  P) 

3n2' 


v  /  32w  32w  ,  0 

K0(°13F  +  °23^  +  2t 


32w  \ 
123?3n; 


Ko(ol  kl  ^2^2  ^ 


(2-22) 
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If  we  define  the  column  matrix  U  as 

U=(J!j  (2-23) 

then,  eqs . (2-21 )  and  (2-22)  can  be  written  in  matrix  form  as 


/l  Ox  9U  .  /0  0\  ||  _  /  0  lx  /32U  32U \ 
O  1  3T  +  (0  c>  u  ■  (-l  0  '3fT  *  17P 


u  w  ,  u  W  ,  o 

—  *  °2W  +  2t> 


0 


92 


w 


°lkl 


a2k2  + 


(2-24) 


Formulation  of  the  Finite-Difference  Equation 

If  we  restrict  ourselves  to  the  first  mode  type  response  of  shells, 
it  is  only  necessary  to  consider  one-quarter  of  the  shell  because  of  the 
symmetry.  Let  us  consider  the  rectangular  network  as  shown  in  Figure  2-3 
at  time  x^.  The  grid  dimensions  in  the  C  and  n  directions  and  the  time 
increment  are  denoted  by  AC>  An  and  Ax,  respectively.  If  one-quarter  of 
the  shell  is  divided  into  (M-l)x(N-l)  sub-divisions,  the  spacing  dimen¬ 
sions  A^  and  An  are 


AC  = 


X 

2 (M-l ) 


An 


1 

2TNT77 


(2-25) 


y  *  A?/An  =  X  (N-1 )/ (M-l ) 

Besides  the  interior  points  (M-l)x(N-l)  and  the  boundary  points  (M+N-l), 
fictitious  points  are  introduced  along  the  lines  C  =  -  AC  and  n=  -  An.  Let 
us  denote  any  variable  ¥  at  a  discrete  point  P.  .  and  at  time  level  x.  to 

I  9  J  K 

k 

be  f.  ..  Hereafter,  the  subscript  represents  the  position  and  the  super- 
script  denotes  time  level. 
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We  shall  reduce  eq - (2-24 )  to  the  finite-difference  form  using  Crank- 
Nicolson  finite-difference  scheme.  The  partial  derivatives  are  approximated 
by 

/3U\  _  /. ,k+l  * | k  \/A 

(57}1  .j  •  (ui ,j  '  ui.j)/iT 


/9U\  _  /..k+1  nk+l  nk  Mk  \/AAcr 

(3T>1.j  -  (u1+u  -  u1-u  +  u1+l,j  -  ui-lj>/4A« 


\  -  (nk+l  _2u^+^  +  uk+1  +  Uk  -2Uk 

W'U  u1+l,j  2Ui  ,j 


+  i -1 »j )/^(A^)2 


(2-26) 


(U).  .  =  (Uk+^  +  Uk  ,)/2 

1  >J  1  >J  >  io 

k  k 

and  (w).  .  can  be  expressed  in  terms  of  vq  .  and  1C  .  as  follows: 

‘«>u-  + 

v?.j  ■  <J  - 

From  these  two  equations 


Using  eqs.(2-26)  and  (2-27),  eq.(2-24)  is  approximately  written  as 

,0  0  x  /..k+1  .  ,,k  x_/  0  1\  1 


Uk+!  -  Uk 


/I  Ox  i  ,.i  wi,j  .  ,0  0  x  (uk+l  x  yk  x  /  u  lx  I _ 

(0  1}  At  (0  c/2}  {Ui,j  Ui , j  l-l 


(2-27) 


+  <UU+1  +  “ij-l 


.k+1 


k+1 


.k+1 


.k+1 


n.ik  x "]  -  /0  Ox  o  i  ,  j  /i, k+1  ,  ,|k+l  9i|k+l  ^ 

“  2Ui,j)]  +  (0  i)-7PnT?p(u1+l,j  Ui-1  ,j  2Ui> y 


,  (0  (Uk+1  +  ||I<+'  .  2UW)  +  (0 

+  (0  1 1  2  ifli lu1  ,j+l  ui,j-l  i  »j  *0  lMTAnl^'V 
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(2-28) 
cont 1 d 


x  (Uk+1  -  Uk+1  -  Uk+1 

x  ui+l  ,j+l  Ui-1  »j+l  ui+l  ,j-l 


+  Uk+]  .  .)  +  (°) 

i-l.j-1'  VA; 


where  k 

A  =  yo.f^k,  +  ®2i  ,jk2  +  o  +  — -f4(w 


(An)2r2  1+1 J 


02 


+  wk 


-  2wlj>  +  T^W.j+i  +  -  2wi,j> 


1-1  »j 


(2-28) 


+  («i+i  1+, 

2(Sn)2r2  ,J 


k  .  k 


wi-lJ+l  "  wi+l  ,j-l  +  wi-l  ,j-l )] 


(2-29) 


Transposing  the  quantities  at  time  to  the  left  hand  side  and  those 
at  time  to  the  right  hand  side,  eq.  (2-28)  is  written  as 

+  +  +  ^UiJ+l 

+  CD5J  u^j  +  [D4]Uk^.1  +  [D3]Uk:]>j+1  +  CD2] 


+  Ui-1 ,j-l  Ui+1  ,j+l  +  ^C8^  Ui+1 , j 

+  tC7]  Ui+1 ,j-l  +  tC6]  Ui,j+1  +  [C5]  Ui,j  +  tC4]  Ui,j-1 

+  [C3]  Ui-l,j+l  +  [C2^  Ui-l,j  +  Ui-l,j-l  +  ^J?A  ) 


(2-30) 


where  [Dm]  and  [C^]  (m=l ,2, "  * ,9)  are  the  square  matrices  defined  by 
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[Dg] 


Kot12^At^ 

4(An)zr 


(2-31) 


[d8] 


At 


0  -  2(AH)2r2 


Ax 

2Ta^FP 


KQai(An)^ 

2('An)2r^ 


CD,]  = 


W*>‘ 

4(An)zr 


[06] 


At 

2Ta^F 


At 

2(An)z 


"H in)2 


[D5]  ■ 


^  (1+12) 

(in)2  r 


TiTP0*  72> 


)+  Ko°l(M)2  +  Ko°2(At)2  +  it, 


(in)2 


(in)2r2 
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Ax 

ZT^V 


At 

zJmY 


K0a2(Ax)2 

ZWT 


0 


0 


(2-31) 

cont'd 


[d3]  = 


0 


KqT12(At)2 

4(An)^r 


0  “  At 

2(An)zr2 


[d2] 


At 

J('An)V 


Kqo-|  (At)2 
2 (An )  zr2 


fT 

0 


0 


CD13 


0 


Kot12(at)2 

4(An)2r 


0 


At _ 

2(An)2rz 


[C8]  = 


At 

2 ( At )zrz 


0 


19 


[c5] 


Tf^F(,+P_) 


iT  (i+  4) 


1- 


AtC 


[C4]  = 


—At 

JJmT2 


At 
2 ( An 


(2-32) 


[c6]  -  [c8] 

[c2]  =  [CQ] 

[Cg],  [Cy] ,  [C3]  and  [C-j ]  are  null  matrices. 

The  compatibility  equation  (2-9)  at  time  t.  and  at  interior  point  P.  . 

K  1  9  J 

(3  <^i  <_M+1 ,  3  <_j<_  N+l )  has  the  following  finite-difference  approximation: 
2Fi  ,j<3r2+4+  F^)+2Cf1+1  >jtl+Fl-l  ,j+l+F1-l  ,j-l+F»l  ,M> 

*  r^Fi+2 , j+l"2Fi  ,j+2+  72Fi-2,j+r2Fi,j-2-4(1+  P^i+l  ,j 


•4(,+r2)Fi  ,j+r4(1+  p)Fi-i.r4(r2+1)Fi.j-i 

=  TC^I+l  ,j+rw1-l  ,j+l+wi-l  ,j-1"wi+l  ,j-l )2 


-<Vi  J*1-1  ,r2wi .j’  (wi .j-r2wi .j> 
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(2-33) 

cont'd 


•ki(wi,j+i+wi,j-r2wi,j)  (4n)2  •k2<wi-n,j+wi-i.j'2wi,j)  (4n)2 


Now  let  us  formulate  the  boundary  conditions  for  the  simply  supported 
shells.  For  the  point  P..  ^  and  P2  m  on  the  boundary,  eqs.  (2-10)  and  (2-11) 


are  written  in  the  finite  difference  form  as 


wi  ,2 

=  0  and  w.  9  =-w .  , 

1  ,0  1,1 

2  £  i  <  M+l 

at  n=0 

Wo 

=  0  and  w9  m  =-w-,  m 
3,m  1  ,m 

2  £  m  <_  N+l 

at  5=0 

2  ,m 

r2(Fi 

+F,  -2F  2)-v(F,+, 

,2tFi-l,2-2Fi 

,2>=0 

(2-34) 


2  <  i  <  M+l  at  n=0 


(F,  +F-,  -2F9  )  -vr2(F9  ,  +  F9  ,  -  2F9  )  =  0 

3,m  1  ,m  2, nr  2, m+l  2,m-l  2, nr 


2<m<N+l  at  5  =0  (2-35) 

The  integrals  in  eqs.  (2-13)  and. (2-14)  are  approximated  by  the  trapezoidal 
rule  as  follows.  In  the  £  direction,  keeping  n  constant,  eq . (2-13)  is 
approximately  expressed  by 

M-l 


a2F  a2F, 


3  £  j  £  N+l 

Similarly,  in  the  n  direction,  eq.  (2-14)  is  expressed  by 


N-l 


5T— T  \  / 9 2F  92F  \  ,  ,9"|-  9"h  \ 

^  i(3F‘  Vij  +  {W  'va^i,j+i 

j=2 


f  9  2F  92F 1 


(2-36) 


3  £  i  ^  M+l  (2-37) 

Equations  (2-36)  and  (2-37)  are  expressed  by  the  following  finite-difference 


equations : 


M-2 

i=2  l  1s 


j+1  +  Fi »j-l 


-2Fij)-^Fi+lJ  +  Fi-U 


"  2Fi,j)  +  r  (F1+i  j+i  +  Fi+i,j-i  '  2Fi+l,j) 


"w1-l,j^  "  (An^klwi,j  +  4(wi+2,j  “  wi,j^ 


M-2 

-  -(Fi+2,j  +  Fu  -  2Fi+i,j>}  =27 

i=2 


-  (in)  knw1+) ,j] 


3  <  j  <  N+l 


N-2  ,  , 

^^<Fi 

j=2 


i+1  ,j  +  Fi-1  ,j 


-  2Fi,j>  -  -<Fi,j+l  +  Fi,j+1> 


-  +  +  Fi-i.J+i  ■  2Fi ,j+i > 

-“<Fi,j+2tFi,j-2Fi,j+l>} 

N-2 

[l(wi,j+l  ‘  wiJ-l)2  ‘  (An)2k2wi,j 


+  ?(wi,j+2  '  wi,j)2  -  (An)2k2wi,j+l] 


(2-38) 


3  <  i  <  M+l 


(2-39) 


22 


For  the  clamped  shells,  the  finite-difference  forms  of  eq .  ( 2-1 5 )  are 


1.2  -° 

andwi,3  =wi,l 

2  ^  i  £  M+l 

at  n  =  0 

(2-40) 

2  ,j  =  0 

andW3,j  -  W1.J 

2  £  j  £  N+l 

at  £  =  0 

(2-41 ) 

Equations  (2-38)  and  (2-39)  must  hold  for  this  case,  too. 

Procedure  for  Solving  w.  .  and  F.  . 

 1  >J  1 

Together  with  eqs.(2-33),  (2-38)  and  (2-39),  the  following  simultaneous 
equations  for  the  stress  function  F  at  the  interior  points  and  the  boundary 
points  are  formulated: 

[A]  {F}  =  {G>  (2-42) 

where  [A]  is  the  (M+N+MN-1 )x(M+N+MN-l )  square  matrix  whose  elements  are 

constants,  { F}  is  the  column  matrix  of  F.  .  and  {G }  is  also  the  column  matrix 

■ 

whose  elements  are  functionsof  W.  ..  In  the  above  formulation,  the  stress 

■ *  *J 

function  at  the  corner  F^  ^  bas  been  assumecl  to  be  zer0  following  Wang's 
assumption  [26].  From  eq.(2-42),  if  the  deflection  w.  .  at  each  discrete  point 
is  given,  then  the  stress  function  F.  .  and  the  membrane  stresses  at  each 
point  can  be  determined. 

The  initial  conditions  (2-26)  now  take  the  form 

w?  .  =  0  and  U?  .  =  0  (2-43) 

1 

Equation  (2-30)  together  with  the  boundary  conditions  (2-34)  or  (2-41)  may 

★ 

be  solved  for  U.  .  using  a  digital  computer  by  the  following  procedure. 

Step.l  Given  the  initial  conditions  (2-43),  eq.  (2-42)  is  solved  for  F?  • 

•  »J 

from  which  the  membrane  stresses  at  each  point  P.  •  are  determined  as 

■  >J 
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°li >j  = 
°2i  ,j  = 

T12i,f 


(Fi  ,j+l+Fi,j-r2Fi,j)/(4n)2 

(FHU+Fi-l,j-2Fi,j)/r2<^2 
^Fi+l,j+l  Fi-1  ,j+l "F1+l ,j-l 

+F1-1,J-1)/4(An)2r 


(2-44) 


(2-45) 


Step. 2  Substituting  these  initial  deflections,  the  membrane  stresses  and  Q 
into  eqs.(2-29),  (2-31)  and  (2-32),  eq.(2-30)  represents  the  simul¬ 
taneous  equations  for  U.  .  Solving  for  U.  .  w.  .  is  then  found 

1  >J  •  i  »J  »  i  »j 

from  the  relation 

W',j=Wi,j  +iT(0,)U1.j  <2-46> 

Step. 3  Substitute  w^.  j  into  the  right  hand  side  of  eq.  (2-42)  and  solve 

f°r  ,*•  Calculate  the  membrane  stresses  by  eq.(2-43).  Return 
to  step  2  and  continue  the  procedure. 


2.1.2  Examples 

As  an  example  of  the  present  analysis,  the  steady  state  response  of 
the  following  structures  to  various  magnitudes  of  step  load  were  determined. 

1.  Simply  Supported  Rectangular  Plate  (Aspect  ratio  A=1 ,  1.35,  1.5 
and  2.0) 

2.  Clamped  Rectangular  Plate  (Aspect  ratio  A=1 ,  1.25,  1.5  anc  2.0) 

3.  Simply  Supported  Double  Curvature  Shell  (Aspect  ratio  A=1 ) 

4.  Clamped  Double  Curvature  Shell  (Aspect  ratio  A=1 )  As  previously 
mentioned,  the  supports  do  not  permit  in-plane  displacement  of  the  edges  of 

the  structure.  For  a  rectangular  plate  with  aspect  ratio  A,  the  spacing  dimension 
AC  for  fixing  N=4  and  M=5  was  determined  by  (See  Table  2-1) 


N-l 
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=  A/8  (2-47) 

If  we  choose  M  and  A 5  such  that 

M  =  [1+  X  (N-l)]  (2-48) 

.  £  _  X 

’  2TJTTT  <2-49> 

the  solutions  will  be  more  accurate  than  those  obtained  by  eq.(2-47).  However 
a  much  longer  computing  time  is  required.  Therefore,  the  choice  of  AC  was  made 
using  eq . (2-47 )  by  sacrificing  some  accuracy  of  the  solutions. 

Since  there  is  no  criterion  for  determining  the  time  increment  At  to 
bring  about  numerical  stability  and  convergence  of  the  solution.  At  was 
determined  such  that  further  deceasing  At  did  not  apprecially  affect  the 
response.  For  the  following  values  of  At,  the  numerical  solutions  were 
apparently  convergent: 


Simply  supported  square  plate 

At<5x10  ^ 

Clamped  square  plate 

At£5x10'3 

Simply  supported  rectangular 
plate  with  aspect  ratio  x=2 

At <7x1  O'3 

Clamped  rectangular  plate 
with  aspect  ratio  x=2 

Atj<7x10"3 

For  all  cases,  it  was  found  that  At=0.002  is  sufficiently  small.  Through  the 
present  Chapter  At=0.002  was  used. 

In  Table  2-2  the  results  are  listed  together  with  Kornishin's  static  solutions 
[27]  and  the  difference  e  between  the  results  obtained  in  the  present  study  and 
his  results  are  also  listed  for  comparison.  It  can  be  seen  that  even  if  the 
interior  points  are  9(N=4  and  M=4)  or  12  (N=4  and  M=5)  in  number,  reasonably 
good  results  were  obtained.  In  Figure  2-5,  a  typical  deflection  response  curve 
for  the  case  of  a  clamped  square  plate  is  shown. 
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TABLE  2-1 


CHOICE 

OF  r, 

M  AND  N 

X 

M 

N 

r=A5/An 

1 

4 

4 

1 

1 .25 

5 

4 

15/16 

1 .5 

5 

4 

9/8 

2 

5 

4 

3/2 

Table  2-2 i  Comparison  of  Steady  State  Response  of  Present  Analysis  with 
Static  Analysis  due  to  Kornishin 
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2.2  Circular  Plates  and  Shallow  Spherical  Shells 

2.2.1  Analysis 
Governing  Equations 

The  finite  amplitude  vibration  of  a  thin  elastic,  shallow  spherical 
shell  of  thickness  h.  Young's  modulus  E,  Poisson's  ratio  v,  constant  curva¬ 


★ 


ture  k  and  base  radius  R  after  assuming  axi symmetry  is  governed  by  the 
following  two  simultaneous  partial  differential  equations  [25]  in  the  polar 
coordinates  system  shown  in  Fig.  2-2. 


h92w  4.  u  *  .  n  4  * 

ph— ^2  +  phc  +  Dv  w 


+  k*)  +  q  *  (r,t)  (2-50) 


(2-51) 


where 


v4 


* 


w 


c 


=  the  normal  displacment  of  a  point  in  the  middle  surface 
=  viscous  damping 


F 


=  stress  function 


★ 

q  (r,t)  =  lateral  load 


D 


=  Eh3/! 2(1 -v2) 


The  membrane  stresses  in  the  middle  surface  are  expressed  by 


* 


(2-52) 
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The  radial  strain  e  and  the  circumferential  strain  eQ  are 


i JL.  .  k*w*  + 

9r  K  w  2l9r  ' 


£ 


e 


★ 

u 

r 


★  * 
k  w 


(2-53) 


where  u*  is  the  radial  displacement  of  a  point  in  the  middle  surface  of  the 
shell . 

The  membrane  stresses  are  also  expressed  in  terms  of  er  and  eQ  as 


1- 


V(er  +  ve0} 


(2-54) 


°e  =,-V(ee  +ver) 

1  “V 


From  eqs.  (2-52),  (2-53)  and  (2-54) 


?  * 

u*  , *  *  ,  1 /&  F 

f  ■ k  w  +  e<^7 


9F  , 
Jr9r' 


(2-55) 


*  *  * 

9u  , *  *  1 /  9w  \2  J  1  / aF  3  F  s 

—  =  k  w  -  7(— )  *  £<757  -  v  -r?) 


Bending  moments  are  given  by 
h!  =  -  D(^C  ♦  7  §£> 


3r 

★ 


M*  =  -  D(9W 


2  * 

+  1^) 


(2-56) 


e  w'r9r  '  9r2 


In  the  present  study,  the  following  two  edge  conditions  will  be  considered. 

(a)  Simply  supported  immovable  edge 

(b)  Clamped  immovable  edge 

where  again  the  edges  cannot  approach  one  another.  Let  us  formulate  these 
boundary  conditions  analytically. 
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(a)  Simply  supported  shallow  spherical  shell  with  immovable  edge 

★  ★ 

Since  deflections  w  and  the  bending  moment  are  zero  at  the  edge, 


from  eq . ( 2-56 )  we  have 


2  *  * 
*  n  .  3  W  ,  V  3W 

w  =  0  and  — %  +  -  — 


0  at  r  =  R 


(2-57) 


Also  the  radial  displacement  u  at  the  edge  is  zero. 
From  eq. (2-55) , 

i-  kV  ♦£<!£-  *|£>-0.tr.R 


r3r; 


But,  since  w  =  0  at  r=R, 


at  r  =  r 


(2-58) 


Considering  eq.(2-52),  an  additional  restriction  for  F  is  that  at  the  center 

* 

of  the  shell,  the  membrane  stress  must  be  finite.  Thus, 

3F  -  0  at  r  =  0  (2-59) 

*  9r 

Since  u  =  0  at  r=0  and  r=R,  integration  of  the  second  equation  of  (2-55)  from 
r=0  to  r=R  leads  to 


?<!r>2  +  r<?5F-  = 0  (2-60) 

which  will  be  used  for  the  formulation  of  the  finite-difference  equation. 

(b)  Clamped  shallow  spherical  shell  with  immovable  edge 
The  deflection  and  the  slope  at  the  edge  are  zero. 


Thus 


w*  =  0  and 


3W 

3r 


=  0  at  r=R 


(2-61) 


Conditions  (2-58),  (2-59)  and  (2-60)  must  hold  for  this  case  too. 

For  convenience,  let  us  nondimensional ize  all  equations  by  introducing 


the  following  dimensionless  parameters: 


*  4  4 

q  =  q  R  VEn 
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w  =  w  /h 


s  =  r/R 
F  =  F*/EIV 


c  = 


4  V  ? 
t(D/phR j 

* ,  4  ,  ^ 

c  (PhR4/D) 


k  R2/h 


*  2  2 
o=o  R  / En 
s  r 

*  2  2 
°e  /Eh 


Ko=  12(l-v  ) 


Using  eq.(2-62),  eqs.(2-50)  and  (2-51)  become 

9W  ,  „4  v  r  9 2  F  / 1  ,  9w  \ 

+  C9t  V  W  =  K"^(k  + 


9t“ 


0L9S' 


S  9S  ' 


92W 


92W  _9W 

9S 


+if<k  +  3S 


•  92W 


+  q] 


.  k(£-3  +  £SL) 

S  9S  '9S2  S9s' 


respectively.  The  boundary  conditions  (2-57)  and  (2-58)  for  simply 
edges  are  now  expressed  by 

at  s=l 


n  j  92w  .  v  9w  _  n 
w  =  0  and  3-r  *  j  35  -  0 


92F  9F 


=  0 


9SZ  S9S 

respectively.  Conditions  (2-59)  becomes 

Equation  (2-60)  is  transformed  to 


S'. 


<ffs  -  4?^  i  rl#2  - 


kw]  ds 


at  s=l 


at  s=0 


0  OJ  ^0 
The  boundary  conditions  for  the  clamped  edge  (2-61)  become 


(2-62) 

(2-63) 

(2-64) 

supported 

(2-65) 

(2-66) 

(2-67) 


(2-68) 
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w  =  0  and  ||  =  0  at  s=l  (2-69) 

The  initial  conditions  imposed  here  are 

w  =  0  and  p-  =  0  at  t=0  (2-70) 

The  problem  for  a  simply  supported  edge  consists  in  finding  w  and  F  from 
eqs.(2-63)  and  (2-64)  subject  to  the  boundary  conditions  (2-65)  through  (2-68) 
and  initial  conditions  (2-70).  For  a  clamped  edge  the  boundary  conditions  are 
given  by  eq.  (2-69)  instead  of  eq . ( 2-65 ) .  The  method  of  solution  employed 
here  is  the  same  as  that  used  for  the  rectangular  shell. 

To  reduce  the  order  of  the  derivatives  in  eq.  (2-63),  we  introduce  two 
new  variable  W  and  V  defined  by 


(2-71) 


V  =  is 
v  9x 

Eliminating  w  from  above  two  equations,  we  have 

9W  _92V 
9t  “9S7 

Using  W  and  V,  eq.(2-63)  is  expressed  by 


w  +  <*  ■  -< 


92W 

FT2’ 


2  9W  _  W  ,  1  9Wx 
?  9?  P+  ?39?} 


(2-72) 


+  K0Cas(W  +  1)  +  oB(^+  k)  +  q] 
If  we  use  the  vector  U  defined  by 


eqs.(2-72)  and  (2-73)  are  expressed  in  the  following  matrix  form: 


'1 

o' 

3U  + 

0 

o' 

u  =  - 

'0 

*1 

1 

92U  , 

'0 

0' 

.0 

1 

F  + 

.0 

.1 

0 

2 

-s 

0 

/ - 

o 

o 

\ _ 

-«■  >•» 

0 

U  + 

/oq+(l Sv  +  V<Y%  >k  , 

(2-73) 

(2-73) 


(2-75) 


Formulation  of  the  Finite-Difference  Equations 
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Consider  a  system  of  one  dimensional  equally  spaced  discrete  points 

k 

P.j  at  time  x^  as  shown  in  Figure  2-4.  Here  the  subscript  i  represents 
the  spatial  position  and  the  superscript  k  denotes  the  time  level.  Be¬ 
cause  of  the  symmetry  we  consider  the  motion  of  only  a  radius  of  the 

shell.  Let  us  divide  the  radius  1  into  N  intervals  and  call  the  point 

k  k 

at  the  center  of  the  shell  P  .  One  fictitious  point  PN+1  outside  the 
edge  is  introduced.  The  spacing  dimension  as  is  1/N. 

Let  us  write  the  finite-difference  analog  for  eq . ( 2- 75 )  using 
the  Crank-Nicolson  finite-difference  scheme  which  was  introduced  in 
the  previous  section.  Since  we  have  a  singularity  at  the  center  s=0, 
we  must  consider  the  case  for  s=0  spearately. 

|y 

At  points  P.j  (0£i<N),  eq .  ( 2-75)  may  be  approximated  by  the  following 
finite-difference  expression: 


tatKok<°0i+‘si> +  (-TZI 


w^  _wk  (Vk+1-Vk+1) 

wi+l  1-1  +  At  Vi+1  i-l; 


4as 


(2- 


l  i 

Here  it  has  been  assumed  that  the  membrane  stresses  remain  constant  in  th 
time  interval  [t^^^+AtJ.  Rewriting  eq.  (2-76),  we  have 
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[b3]<;]  +  [B2JU*5+1  +  [B,  ]u^| 


■  ^“l+l  +  ^  *  ^-l  + 


where  [B.]  and  [A.]  (j=l,2,3)  are  square  metrices  defined  by 

J  J 


(2-77) 


[b3] 


0 


At 

'2T4TP 


At  (At  )2  /I  i/  k  \ 
2ass.  4ass.  's?~oCT0i' 

1 


[b2]= 


At 

(as)2 


to VsV  h  +  V&  1+I-J 


(2-78) 


[B-,] 


0 


At 

'  2Ia?F 


A(1  .1) 

2asvas  s..' 


(Atf  fl 
4ASS.j  VSi 


-K 


k  \ 
oa0i 


0 


[A3] 


At  /I 
2AS^AS 


+ 


At 

2(as )z 


0 


1 


At 

T 


(K  acK,+ 
v  o  si 


tttt 


s) 


At 

TFT2 

1  CAT 
2 


[A,] 


0 

At  /  1 

2as's^ 


At 

FTasT2 


0 
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k 

and  G.  is  defined  by 

<  •  K0Axtq+K(oekjW^)]  ♦  *x(°|l  - 

,k  ; 


The  compatibility  equation  (2-64)  at  the  point  (l£i<N-l)  may  be 
written  in  the  finite-difference  form 

(Fi+2”4Fi+l  +6Fi”4Fi  -1+Fi -2>/^s>2  +  (Fi+2-2Fi+l+2Fi-l 

-Fi_2)/Si(As)3-  (Fi+1-2Fi+Fi.1)/(siAs)^  +  (F.^-F..-,) 

x(l/2Ass|)  =  -  (wi+1-2w.+w._1)(wi+1-wi_1)/2si(As)3 

-  k(wi+1-2wi+w._1 )/ ( AS )2  -k(wi+1-wi_1)/2si(As) 

At  the  point  Pk,i.e.,  the  center  of  the  shell,  using  L1 Hospital's  rule, 
eq.(2-75)  takes  the  form 


(2-79) 


f 

1  0 

M  + 

f 

0 

0 

u  =  - 

0 

-1  ' 

92U 

+  K0(a$+ae) 

f 

0 

o' 

9t 

9SZ 

0  1 

0 

c 

8 

0 

1 

0 

**  / 

y 

Lj 

9 

\ 

•Ko{q+0s+ae) 

which  has  the  following  finite-difference  analog: 
2[B']Uk+1+  [B’]Uk+1  =  2[A^]Uk  +  [A’]Uk 


(2-80) 


Kotq+k(as0+ae0)]AT 


(2-81) 
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where 


[Bj]  = 


At 

2 JJTF 


4A 
.3(as)2 


[BJ] 


4VaeO+CTs(PAT' 


8  At 

3  TasT2 


At 

TaIF 


tjcat  +1 


(2-82) 


At  s=0,  the  compatibility  equation  (2-64)  can  be  expressed  as  (See  APPENDIX), 

(2-83) 


83  F  _  /  32W  qi.32W 

3„  4  "  ’  "  AKisT 


3s 

which  reduces  to 

|(F2-4F1+3F0)  =  -  (wrw0)  -k(wrw0)  (as)2 
The  boundary  conditions  (2-65)  and  (2-66)  are  expressed  by 
w. 


N 


0  and  5!*^N=L  .  0 


(AS T 


sN  2as 


FN+1+FN-T2FN  FN+1"FN-1  n 
-  *"S^s  -  0 


T5F 


(2-84) 


(2-85) 


(2-86) 


respectively.  The  integral  in  eq.(2-68)  is  approximated  by  the  trapezoidal 
rule  as 


1/1  ,,\  /3^.F\  +'4_7/  3F  S2F x  1/3F  32Fv 
2-n-v)  (gp0o  +  >-,  (  S3S  -  3^7) i  +  2^ S3?  ~vdST) 


N 


=  "  kw0]  +  kwN] 


(2-87) 


In  finite-difference  form,  this  becomes 


F,tF0 


*=*■  Fi+rFi-i 


AS 


,  1  (WFN-1 
2  \  2s,,as 


’r 


Fi+i+Fi-r2Fi 

[AiP 

FN+1+FN-1~2FN 

(As)2 
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+ 


i+1 _wi-l 

■  8TasT' 


)2 


1  /l  ^WN+TWN-1^2 
+  2  (l - JJsf~ 


(2-88) 


The  boundary  conditions  (2-69)  for  a  clamped  edge  are  expressed  by 


WN  *  0  and  Vl  *  WN-1  (2-89) 

Also,  conditions  (2-86)  and  (2-88)  must  be  satisfied  for  clamped  shells. 

The  initial  conditions  (2-70)  have  the  form 

w?  =  0  and  V°  =  0  1-0,1 ,••• ,N+1  (2-90) 

Procedure  for  solving  for  and  F. 

Equations  (2-79),  (2-84)  and  (2-89)  constitute  the  following  simultaneous 
equations  for  F.: 

[ K] { F }  =  {L}  (2-91) 

where  [K]  is  an  (N+l)x(N+l)  matrix  whose  elements  are  function  of  as,  and  s^ 
(F }  is  a  column  matrix  of  F^ ,  and  {L }  is  a  column  matrix  whose  elements  are 

k 

function  of  w^ ,  If  the  deflections  at  every  point  are  known,  then  the 
stress  function  F  at  every  point  can  be  determined  by  solving  eq.(2-91), 
namely 

{F}  =  [K]"1 {L }  (2-92) 

The  procedure  for  solving  for  w^  is  essentially  the  same  as  that  described 
in  the  previous  section  (2-1-1): 

Step  1.  Substitute  the  initial  conditions  (2-90)  into  (L }  of  eq.(2-92)  and 
obtain  F°.  Then  the  membrane  stresses  at  each  point  can  be  cal¬ 
culated  by 
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Fi+rFi-i 

2s  ^  AS 


1  <  i  <  N 


(2-93) 


ei 


Fi+i+Fi-r2Fi 

(AS  J7 


At  the  center  of  the  shell  s=0,  using  L'Hospital's  rule 

°s  =o0  -  2(F1-F0)/(4S)2 
0  0 


(2-94) 


Step  2.  Substitute  a  and  oQ  into  (2-78)  and  (2-82),  and  compute  U?.  Solve 

si  0.  i 

1  1 

eqs.  (2-77)  and  (2-81)  for  U-.  The  deflection  can  be  obtained  from 


w]  =  W°  +  At (0, 1 ) U j 


=  W?  +  At vj 

Step  3.  Substitute  the  new  deflection  wj  into  {L }  of  eq.  (2-92)  and  compute 
f|.  Repeat  this  procedure. 

2.2.2  Examples 

As  an  example,  the  computation  of  the  deflection  response  was  carried  out 
for  the  following  cases. 

1)  Simply  Supported  Circular  Plate 

2)  Clamped  Circular  Plate 

3)  Simply  Supported  Spherical  Shell  (k=l,2,  and  3) 

4)  Clamped  Spherical  Shell  (k-1,2,  and  3) 

where  the  edges  cannot  approach  one  another.  The  spacing  as  was  taken  to  be  0.2. 

The  time  increment  At  for  numerical  stability  and  convergence  of  the  solution  for 

the  given  as=0.2  was  found  to  be 

_2 

At<_1.0x10  for  a  simply  supported  edge 

_3 

At<2.5x10  for  a  clamped  edge 

However,  throughout  this  section  the  deflection  response  was  calculated  using 
At=  0.002.  Typical  deflection  response  curves  for  a  clamped  spherical  shell 
are  shown  in  Fig. 2-6.  The  steady  state  responses  of  the  above  cases  in  the  presence 
of  heavy  damping  are  listed  in  Table  2-3  together  with  Kornishin's  static  solutions. 
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III 

SIMULATION  OF  A  STATIONARY  GAUSSIAN  WHITE  NOISE 

3.1  White  Noise 

White  noise  is  a  mathematical  idealization  of  a  stationary  random 
process  in  which  the  power  spectral  density  is  constant,  GQ  at  all 
frequencies.  Since  the  autocorrelation  function  R(t)  for  white  noise 
is  expressed  by 

R(t)  =  2ttG06(t)  (3-1) 

where  6(t)  is  the  Dirac  delta  function,  the  process  has  an  infinite 
variance  and  is  completely  uncorrelated  at  different  times. 

3.2  Simulation  of  Approximate  White  Noise  Processes 

True  white  noise  is  physically  impossible  and  cannot  be  simulated 

because  it  requires  an  infinite  mean  square  value.  Therefore  in  this 

section,  approximate  Gaussian  white  noise  whose  power  spectral  density  is 

constant  over  the  range  of  frequencies  of  interest  and  falls  to  zero  as 

the  frequency  tends  to  infinity  is  generated.  To  do  this,  first  a  sequence 

of  independent  random  numbers  Bk  distributed  uniformly  in  the  interval 

(0,1)  is  generated  by  a  power  residue  method  [29]  on  the  CDC  3600 

digital  computer.  Then  a  white  sequence  of  Gaussian  numbers  vk  with 

2 

mean  zero  and  variance  ay  is  obtained  through  the  transformations 

1/2 

vk  =  ay(-2  loge  Bk)  cos2ttB 

1/2 

vk+l  ’  °v(-2  lc%  Bk>  s1n2’Bk+l  J 


k  odd 


(3-2) 
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Construct  a  function  p(t)  consisting  of  a  sequence  of  step  functions  having 
constant  time  interval  Ah,  with  the  ordinates  in  the  various  intervals 
being  white  numbers  vk.  The  mathematical  formulation  of  the  function 
p(t)  is 

00 

p(t)  =  Z  v.  { u [ t  -  (k+l)Ah]  -  u(t  -  kAh)}  (3-3) 

k=^o 


where  u(t)  is  a  unit  step  function  defined  by 
u(t)  = 


1  t>0 

1/2  t=0 

0  t<0 


(3-4) 


In  this  study,  the  process  is  assumed  to  be  ergodic,  stationary,  and 

Gaussian.  Hence  ensemble  average  may  be  replaced  by  time  average. 

2 

The  temporal  mean  value  ^  and  variance  an  become,  respectively, 


“p  =  i1m  5T 

r  j  -*» 


1 


-T 

N 


p(t)dt 


(3-5) 


lim^r  z  v, 
N+~  k=-N 

T 

-T 
N 


ap  =  iim  1 T 


T->' 


[p(t)]"dt 


(3-6) 


=  lim 


1 


2  _  2 

w  2N  , Z  Mvk  ar 
N^-“  k=-N 


The  autocorrelation  function  of  the  random  process  p(t)  is  calculated  from 

rl 


R(t)  =  lim 


1 

7T 


p(t)p(t+x)dt 


(3-7) 
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(3-8) 


which  yields 

(  av(1  aTF^-  ^  M  -  Ah 

R(t)  =  < 

|  0  | t  |  -Ah 

The  power  spectral  density  function  G(w)  is  the  Fourier  transform  of 
the  autocorrelation  function  R(x).  Since  R(t)  is  an  even  function  of 
t,  then  we  have 


G( 


u)  =  P  R(x)e  wlT  dt 

j  -00 


r 

j  -« 


R(t)  COSordT 


f'sln^ 

U)Ah 

I 

(3-9) 


For  uAh  small,  G(w)  can  be  approximated  by 

v2 


G(„)  *  l^L  [1  .  is*jjl£] 


(3-10) 


Choosing  Ah  sufficiently  small  G(w)  can  be  approximated  by  o^Ah/ir. 
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IV 

STATIONARY  RESPONSE  OF  PLATES  AND  SHALLOW 
SHELLS  TO  GAUSSIAN  WHITE  RANDOM  EXCITATION 

4.1  Introduction 

In  the  present  Chapter  the  response  analysis  of  plates  and  shallow 
shells  to  stationary  Gaussian  random  excitation  will  be  carried  out 
by  various  methods.  All  equations  and  nomenclatures  used  here  are  the 
same  as  those  described  in  Chapter  II. 

The  random  excitation  q(£,n,x)  or  q(s,x)  is  assumed  to  be  uniformly 
distributed  over  the  structure  and  applied  normal  to  the  middle  surface 
of  the  structure.  Furthermore,  the  random  excitation  is  also  assumed  to 
be  stationary,  ergodic  and  Gaussian  white  noise  with  zero  mean  value,  that  is 

q(5.n,x)  or  q(s ,t)  =  Q(x) 

E[Q(x)]  =  0  (4-1) 

E [Q ( x i ) Q ( X2 ) ]  =  2ttGq6 (ti  -  X2) 

where  GQ  is  constant,  6(t)  is  the  Dirac  delta  function  and  E  [  ]  is  the 
expectation  operator.  The  methods  of  solution  are  the  following. 

Method  I.  Numerical  Simulation  by  the  Runge-Kutta  Method 

The  equation  of  motion  and  the  compatibility  equation  are  reduced 
to  a  single-degree-of-freedom  dynamical  system.  To  do  this  we  assume 
the  normal  deflection  w(c,n,x)  or  w(s,x)  to  be  the  product  of  the  coordinate 
function  ^(c»n)  which  satisfies  the  boundary  conditions  and  is  suitable 
for  representing  the  deflected  shape  of  plates  or  shells  and  the  generalized 
coordinate  f(x)  such  as 
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w(c,n,x)  =  f  (x)ip(c»n) 

w(s,x)  =  f(xMs) 


(4-2) 


Then  we  substitute  eq.  (4-2)  into  the  compatibility  equation  and  solve 
for  the  stress  function  F  in  terms  of  f(x).  Upon  using  Galerkin's 
method  [25],  the  equation  of  motion  and  compatibility  equation  yield 


f  +  2cuof  +  *2(f  +  ef2  +  Yf3)  =  <xQ(t)  (4-3) 

2 

where  c  is  the  fraction  of  critical  damping.  u)Q,  e,  y  and  a  are 
constants.  Note  that  for  a  flat  plate  6  is  zero. 

The  approximate  white  noise  process  p(x)  is  generated  digitally 
by  following  the  method  described  in  Chapter  III,  by  use  of  eqs. 

(3-2)  and  (3-3).  The  power  spectral  density  is  expressed  approximately  by 

G(«0*  ^  [1  -  %^2]  (4-4) 


If  Ah  is  small  and  the  damping  of  the  system  x,  is  small,  G(u>)  is  almost 
flat  near  the  frequency  coQ  and  p(x)  can  be  considered  to  be  white  noise. 
Therefore, 

G(«)  *  a2  Ah/ir  =  G0  (4-5) 

The  magnitude  of  the  power  spectral  density  can  be  varied  by  changing 
either  ay  or  Ah.  Equation  (4-3)  is  now  integrated  numerically  by  the 
Runge-Kutta  method  employing  the  p(x)  from  eq.  (3-4)  as  an  input  Q(x). 

In  the  present  study,  the  integration  time  increment  Ax  and  the  basic 
time  increment  Ah  are  chosen  as 


Ax  =  Ah  =  1/1  OgOq 


(4-6) 


47 


Since  the  random  process  is  assumed  to  be  ergodic,  the  ensemble 
average  may  be  replaced  by  the  temporal  average.  Therefore  the  mean, 
mean-square,  and  the  variance  of  the  response  f(x)  are  computed  by 


E[f]  * 


1 

T 


rT 

f  dx 

J0 


1 

N 


N 

E 

k=0 


f 


k 


E[f2]* 


1 

N 


N 

E 

k=0 


f 


2 

k 


2 

af 


1 

N 


N 

E 

k=0 


N 

E 

k=0 


fk> 


(4-7) 


2 

where  is  the  variance  of  f(x)  and  f^  is  the  response  of  f  at  time  x^ 

2  2 

The  mean-square,  E[w  ],  and  the  variance,  a^,  of  the  central  deflection 
are  found  from  eqs.  (4-2)  and  (4-7)  as 


E[w2]  =  E[f2]  /  (1/2 X,  1/2) 
°w  =  °f  ^  0/2  A,  ]/2) 


(4-8) 


The  averaging  time  T  (or  N)  in  eq.  (4-7)  must  be  determined  such  that 

O 

further  increasing  T  does  not  affect  E[f]  or  E[f  ]  appreciably.  Here  T 
is  chosen  to  be  6600  Ah. 

Method  II.  Numerical  Simulation  by  the  Finite-Difference  Technique 

Insteady  of  solving  the  reduced  single-degree-of-freedom  dynamical 

h 

equation  (4-3),  the  equation  of  motion  (2-8)  and  the  compatibility  equation 
(2-9)  are  integrated  numerically  by  the  finite-difference  method  employing 
the  simulated  white  noise  p(x)  as  an  input.  The  integration  time  increment  Ax 
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and  the  averaging  time  T  are  chosen  to  be  the  same  as  those  used  in  method 
I.  The  mean-square  response  of  the  deflection  is  then  calculated  by  use 
of  eqs.  (4-7)  and  (4-8) . 

Method  III.  Use  of  the  Fokker-Planck  Equation  [14] 

Since  the  input  Q(tO  in  eq.  (4-3)  is  white  noise,  this  equation  can 
be  solved  exactly  by  the  Fokker-Planck  approach.  Writing  y-j  =  f  and 
y2  =  f,  eq.  (4-3)  is  equivalent  to  the  following  pair  of  first-order  equations. 

*1  =  y2  (4-9) 

o  o  o 

y2  =  -2;u>0y2  -  wQ(y-|  +  ey-j  +  yy-j)  +  Q(x)a 


The  stationary  Fokker-Planck  equation  associated  with  eq.  (4-9)  is 
2 

~T~  Z2  "  l^(y2p)  +  dh  {[2?V2  +  wo(yl  +  6yl  +  y?y)]}  p 


ay; 


ay2 


=  o 


(4-10) 


where  p  is  a  joint  probability  density  for  y-j  and  y2«  The  solution  for 
p(y1 ,  y2)  obtained  by  Caughey  [14]  is 


p(yr  y2)  =  p(f»  f) 


cexp  {-  []■  f2  +  1  ef3  +  \  yf4  +  ^)2]}  (4-11) 

ira  Gq  wO 


where  C  is  a  normalization  factor  defined  by 


C  =  l/[ 


p(f ,f )dfdf] 


(4-12) 


Since  f  and  f  are  statistically  independent. 


P(f.f)  =  P(f)p(f) 
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where  p(f)  is  expressed  by 


P(f) 


exp  ( -c-| f*"  -  c2f3  -  Cgf^) 

exp  (-c-|f^  -  c2f^  -  c-jf^df 

—00 


where 

c}  =  1/(2 o2f  ) 

0 

C2  =  6/ (4a?  ) 

0 

c3  =  y/(4a?  ) 

J  To 

°f  =  7tG0°‘2/(44Uq) 


2 

and  af  is  the  mean-square  response  of  the  linear  system 
0  2 

to  eq.  (4-3).  The  variance  of  of  f  is  obtained  by 
4  =  E[f2]  -  (E[f])2 

■  f2p(f)df  -  [ r fp(f)df]2 

J  .00  J.  00 

p 

For  plates,  i.e.,  e=0  and  E[f]=0,  af  can  be  expressed  by 
cylinder  function  D^(z)  as  [16,28] 

4  =  E[f2] 
r  2 

=  rp(f)df 

J  -<x>  I top  r,  . 

_af^’exp(-c-| t  -  c3r)df 
°°exp(-c-|f^  -  c3f  )df 


(4-13) 

(4-14) 

=y=0)  corresponding 

(4-15) 

parabolic 


50 


Introducing  the  change  of  variable  =  g,  we  find: 


g1/2exp(-c]g  -  c3g2)dg 


exp( -c-, g  -  c3g  )dg 


J3/2 


^2c3^ 


ra 


2(2c3) 


w 


1/2 


(2c,) 


W 


(4-16) 


In  the  present  study  the  integrals  of  eq.  (4-15)  are  evaluated  numerically. 
Method  IV.  Use  of  the  Equivalent  Linearization  Technique. 

Let  us  assume  that  the  nonlinear  differential  equation  (4-3)  is 
linearized  as 

f  +  2c«0f  +  ^f  =  «Q(t)  (4-17) 

where  m  is  the  equivalent  linear  stiffness.  The  error  e  caused  by 
oe 

this  linearization  is  the  difference  between  eq.  (4-3)  and  eq.(4-17),  i.e., 

e  ■  (“oe  -  "o)f  '  “o(ef2  +  Yf3)  <4'18) 

2  2 

In  order  to  determine  uQe>  we  choose  u)Qe  so  as  to  minimize  the  mean-square 

value  of  the  error  e.  The  mean-square  of  e  is 

r-r  2-i  /  2  2\2crjr2-|  ,  4Q2rr.f4-i  ,  4  2r-r.p6-| 

E[e  J  =  (toQe  -  w0)  E[f  ]  +  ioQe  E[f  J  +  wqy  E[f  J 

-  2(%e  “  wo^o6E^f3-^  +  2woY6E^f5-^  '  2t0oY^woe  ~ 

p 

Minimization  of  E[e  ]  is  achieved  by  requiring  that 

d(“oe> 


(4-19) 
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from  which  we  have 


0)  =  CO 

oe  o 


fl  ,  BE[f3]+yE[f4] 
E[t  ] 


■} 


(4-20) 


Since  Q(t)  is  a  stationary  Gaussian  random  process  with  zero  mean 
value,  the  response  of  the  linearized  system  is  also  assumed  to  be 
stationary  Gaussian  if  the  nonlinearity  of  the  system  is  small.  Then, 
E[f4]  and  E[f^]  are  evaluated  as 


E[f4]  =  3(E[f2])2 


E[n  =  0 


(4-21) 


The  mean-square  response  of  the  linearized  system  is  found  from 


E[f2]  ■ 


2r 
CL  b 

0 

71  2\2 . .  2  2  2 

Ke'V  +4C  V 


dco 


2r 
ira  G. 


2 

4cto.aj„^ 

o  oe 


Substituting  eqs.  (4-21)  and  (4-22)  into  eq.  (4-20),  we  obtain 


(4-22) 


3v(E[f2]2  +  E[f2]  -  4o  "  0 


(4-23) 


where  is  the  mean-square  response  of  the  corresponding  linear  system 

2  2  2 

of  eq.  (4-3)  and  obtained  by  replacing  uQe  by  uQ.  Solving  for  E[f  ]  from 
eq.  (4-23) ,  we  find 


1/2 


E[f2]  -  ^  [0  +  12v40)  - 


(4-24) 


52 


3 

Note  that  if  6  is  not  zero,  E[f  ]  is  not  zero.  However,  the  estimation 
of  E[f  ]  from  the  linearized  system  leads  to  the  conclusion  that  E[f  ] 
is  zero.  Therefore,  for  the  case  of  shells  where  6  is  non-zero,  this 
method  will  not  give  a  good  approximate  solution. 

4.2  Plates 

The  mean-square  deflection  response  is  evaluated  by  Method  I 
through  IV  for  the  following  cases: 

(a)  Simply  Supported  Square  Plate  (x=l) 

(b)  Simply  Supported  Rectangular  Plate  (  x=2) 

(c)  Clamped  Square  Plate  (x=l) 

(d)  Clamped  Rectangular  Plate  (x=2) 

(e)  Simply  Supported  Circular  Plate 

(f)  Clamped  Circular  Plate 

It  is  assumed  that  the  supports  for  all  cases  do  not  permit  in-plane 
displacement  of  the  edges  of  the  structure.  The  governing  equations 


and  the  boundary  conditions  are  described  in  Chapter  II. 

We  assume  the  following  form  for  deflection  w: 

w  =  f‘(x)ip(c,n)  =  f(x)sin^simrn  for  (a)  and  (b)  (4-25) 

w  =  f(T)i^(c»n)  =  f(x)sin2ysin2Trn  for  (c)  and  (d)  (4-26) 

w  =  f(x)iKs)  =  f(x)(l-2m-|S2+m1m2s4)  for  (e)  (4-27) 

w  =  f(x)iHs)  =  f (x) (1-2s2+s^)  for  (f)  (4-28) 


where  m-j  =  (3+v)/(5+v)  and  m2  =  (l+v)/(3+v)  All  above  assumed  solutions 
satisfy  the  boundary  conditions  for  each  case.  By  substituting  w  into  the 
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compatibility  equation  and  solving  for  F,  then  using  Galerkin's  method, 
we  have  the  following  final  equation  for  each  case: 

f  +  2cu)Qf  +  «2(f  +  yf3)  =  aQ(t) 

2  2 
where  u>0>  y  and  a  for  each  case  are  listed  in  Table  4-1.  Once  wQ,  y 

and  a  are  determined,  the  mean-square  response  of  the  deflection  can 

be  evaluated  by  methods  I,  III,  and  IV.  The  damping  coefficient  is 

chosen  as  c  =  0.05  for  each  case.  The  mean-square  response  of  the 

central  deflection  versus  spectral  density  GQ  are  shown  in  Figs.  4-1 

through  4-6. 

The  results  show  that  the  mean-square  response  by  methods  I,  III, 
and  IV  are  reasonably  close  to  each  other,  however  the  results  obtained 
by  method  II  deviate  very  much  from  those  found  by  the  other  three 
methods.  This  is  due  to  the  truncation  error  and  the  propagation  error 
involved  in  the  process  of  numerical  integration.  (Case  (b)  was  not 
worked  by  method  II.)  Since  the  computation  time  required  by  method  II 
is  great  (for  example  for  case  (a),  it  took  46  minutes  to  get  the  single 
point  in  Fig.  4-1.),  the  method  II  is  not  recommended.  It  is  also  to  be 
noted  that  the  results  obtained  by  the  equivalent  linearization  technique 
are  smaller  than  those  by  the  Fokker-Planck  approach  for  all  cases. 

4.3  Shells 

In  this  section,  the  variance  of  deflection  of  the  following  shallow 
shells  is  evaluated  by  methods  I  and  IV. 

(g)  Simply  Supported  Cylindrical  Shells  (k^=0  and  k2=5,  aspect 
ratio  X  =  1.0).  Here,  w  is  assumed  to  take  the  form  of 
eq.  (4-25). 


54 


(h)  Clamped  Cylindrical  Shells  (k-|=0  and  I^B^aspect  ratio 

X  =  1.0).  Here,  w  is  assumed  to  take  the  form  of  eq.  (4-26). 

(i)  Clamped  Cylindrical  Shells  (k-|=0  and  k2=5;  aspect  ratio  x=  2.0). 
Here,  w  is  assumed  to  take  the  form  of  eq.  (4-26). 

(j)  Simply  Supported  Spherical  Shell  (k=0.5).  Here,  w  is  assumed 
to  take  the  form  of  eq.  (4.27). 

(k)  Clamped  Spherical  Shell  (k=1.0).  Here,  w  is  assumed  to  take 
the  form  of  eq.  (4-28). 

Following  the  same  procedure  described  in  section  4.2,  we  have 


f  +  2cu)0f  +  ujj(f  +  Bf2  +  yf3)  =  cxQ(t) 


(4-3) 


where  the  constants  w  ,  B,  y  and  a  for  each  case  are  listed  in  Table  4-1. 

2 

The  variance  o of  the  central  deflection  versus  power  spectral  density 

Gq  for  each  case  is  plotted  in  Fig.  4-7  through  4-11.  From  these  figures 

it  can  be  seen  that  the  result  obtained  by  the  numerical  simulation  method 

2 

are  in  good  agreement  with  the  exact  solution  when  ow  is  small. 
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TABLE  4-1 

VALUES  OF  w3,0  ,  /  AND  a 

1,  simply  Supported  Rectangular  Shells  (Cases  (a)  ,_(b_)_and  ( fQl 


“o 

( x3+l ) 3  T ta/£  +24 ( 1-V3 )  (k1+k2 ) 3  (  tt4-64 )/tt4(  X2+l ) 2 
+768(k3+2Vk1k2X3+lc|x4  )/\4tt4 

p 

|-8(l-v3)[4(k1+k2)/(x3+l)3+(k1+k2x  )/x4]  -72[kx 

+  vx2  (kx+k2)+k2x4]  /x4}/w2 

y 

|3p4[(i-v2)(x4+D/2+(x4+2yx2+D] 

a 

192(l-v2)/n3 

2*  Clamped  Rectangular  Shells  (Cases  (c ) ,  (d )_, ,(h_?  and  Xtll 


“o 

16  [it4(3Xa  +2x3+3)+3(k^+2vx3k1k2+X4k|]/9X4  +(1-V3) 
x  [32(k3+k3X4  )/x4  +  (k1+k2)3/(X2+l)3]  /12 

p 

|-4(l-y3)iT3[(k1+k2x4  )/x4  +  (k1+k2)/(x2+l)2]  -6rr3 
xfk3+yx3(k1+k2)+k2x4]  /x4  J/w2 

y 

{5Ht*(1-v2)(X4+1)/36x,+  |(l-v2)"4[l/(x2+4)2 

+  l/(4x2+l)2  +  4/(x2+l)2]  +  |tt4(x4+2Vx3+D/x*}  A)2 

a 

64 (1-V2 )/3 
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3 •  Simply  Supported  Spherical  Shells  (Cases  (c)  and  ( j  )1 


< 

-  |m2(2m2-3)+y(m|-5m2+5)  (l-v2  )k2-6(2m2-3)  U+V) 

x  [m2(5-V)-3(3-V)]k2]m1A 

3 

5m2-30m2+57m2-36 )k  -  ^  (3m2”^m2"*"^^ 

x [m2(5-V)-3(3-v)]k-  |(l+v)(2m2-3)[m|(7-v)-4m2(5->) 

+  6(3-V)]  kJm2A/aj2 

i 

|~(m2-7m|+266m| -3 5m2+ 14 ) ( 1- y3 ) -  j( 1+ V ) ( 3m|-8m2+  6 ) 
x  [m|(7-y)-4m2(5-v)+6(3-->')]  jrn^A/uj2 

a 

60(1-V3)(7+V)(5+V)/(3V2+36V+113) 

4.  Clamped  Spherical  Shells  (Case  (f)  and _ 

Co2 

0 

4 [320+8  ( 7-2 V)  (l+V)k2] 

3 

-20 (3-V)  (l+l)k/cj3 

r 

”(23-9V)  (1+V  )/co2 

a 

20 ( 1-V2  ) 

where 

A=  -120(3+v) (5+v)/(3v2+36^+113) 
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V 

CONCLUSIONS 

The  large  amplitude  vibrations  of  thin  elastic  plates  and  shallow 
shells  having  various  boundary  conditions  and  subjected  to  random 
excitation  were  investigated  by  using  verious  approximate  techniques. 

As  a  preliminary  study,  the  analysis  of  the  response  of  thin  elastic 
plates  and  shallow  shells  with  damping  undergoing  moderately  large 
deflections  and  subjected  to  step  loading  was  carried  out  by  the  finite- 
difference  method  in  Chapter  II.  The  steady  state  response  for  each 
case  was  compared  with  the  existing  solution  obtained  by  static  analysis 
to  check  the  accuracy  of  the  approximation.  The  integration  time 
increment  bringing  about  numerical  stability  and  convergence  of  the 
solution  for  a  fixed  grid  spacing  was  also  determined. 

In  Chapter  III  a  digital  simulation  technique  for  representing 
a  white  stationary  Gaussian  random  process  was  presented.  The  random 
vibrations  of  thin  elastic  rectangular  plates  and  circular  plates  subject 
to  white  random  excitation  were  simulated  numerically  by  two  different 
methods.  The  first  method  consists  of  three  steps:  The  first  step  requires 
the  reduction  of  the  governing  partial  differential  equations  to  a  single- 
degree-of-freedom  dynamical  equation.  This  was  achieved  by  assuming  the 
normal  deflection  to  be  represented  by  a  product  of  the  generalized 
coordinates  and  a  coordinate  to  simulate  the  random  load  numerically  by 
the  procedure  described  in  Chapter  III.  Then  the  desired  statistical 
properties  were  computed  by  numerically  integrating  the  dynamical  equation 
by  the  Runge-Kutta  method  using  the  simulated  random  load  as  an  input. 
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The  second  method  is  to  integrate  the  governing  equation  of  motion  and 
the  compatibility  equation  numerically  by  the  finite-difference  methods 
described  in  Chapter  II  again  employing  the  simulated  random  load  as 
an  input.  In  the  second  numerical  simulation  method  the  total  response 
may  contain  the  response  associated  with  higher  modes,  however,  this  method 
requires  much  more  computing  time  than  the  first  technique.  To  compare 
the  results  obtained  by  the  two  numerical  simulation  methods,  the  mean- 
square  response  of  deflection  was  determined  by  both  the  equivalent 
linearization  technique  and  by  the  Fokker-Planck  approach.  For  the  cases 
of  rectangular  plates  and  circular  plates,  the  mean-square  response  of 
the  nonlinear  system  was  found  to  depart  more  from  the  linear  response 
for  simply  supported  plates  than  for  clamped  ones.  The  solutions  obtained 
by  the  first  numerical  simulation  method  were  reasonably  close  to  those 
obtained  by  the  equivalent  linearization  technique  and  by  the  Fokker-Planck 
approach.  However,  the  second  numerical  simulation  method  gave  rather 
poor  solutions  because  of  the  truncation  error  and  the  propagation  error 
involved  in  the  integration  process.  It  can  be  concluded  that  considering 
the  accuracy  of  the  solution  and  the  computing  time,  the  second  numerical 
simulation  method  is  not  suitable  for  the  simulation  of  random  vibrations 
of  plates  unless  a  much  more  efficient  method  for  solving  the  governing 
partial  differential  equations  is  used. 

In  the  last  section  of  Chapter  IV,  the  stationary  response  of  shallow 
shells  to  white  random  excitation  was  obtained  by  both  the  first  numerical 
simulation  method  and  the  Fokker-Planck  approach,  within  the  framework 
of  a  single  mode  approximation. 


NOMENCLATURE 


a,  b 

=  dimension  of  plate  or  shell 

[AkL  [Bk] 

=  square  matrices 

[A^],  CBfcl 

=  square  matrices 

[Ck],  [Dk] 

=  square  matrices 

c*,  c 

=  viscous  damping  and  dimensionless  viscous 
damping,  respectively 

D 

=  flexural  rigidity  of  plate  or  shell 

E 

=  Young's  modulus 

F*.  F 

=  Airy  stress  function  and  dimensionless  Airy 
stress  function 

h 

=  thickness  of  plate  or  shell 

{F},  {G} 

=  column  matrices 

{L},  (K> 

=  column  matrices 

kx»  ky»  k* 

k i »  k2 »  k 

=  curvatures  of  shell 

=  dimensionless  curvatures  of  shell 

Ko 

=  12(l-v2) 

Mx,  My 

Mr,  Mg 

=  bending  moments 

=  bending  moments 

M 

=  number  of  interior  points  plus  boundary 
point  in  £  direction 

N 

=  number  of  interior  points  plus  boundary  point 
in  n  direction 

p k  pk 

i,j‘  Pi 

Q 

=  points  at  time  level 

=  amplitude  of  lateral  load 

q*»  q 

=  lateral  load  and  dimensionless  load,  respectively 

r 

=  coordinate 

r 

=  ratio  AC/An 

R 

=  base  radius  of  spherical  shell 

R(t) 

=  autocorrelation  function 

s 

=  dimensionless  coordinate 

AS 

=  spacing  dimension 

t 

=  time 

u(t) 

=  unit  step  function 

U*,  V* 

=  displacements  in  x  and  y  directions,  respectively 

U,  V 

=  dimensionless  displacements  in  c  and  n  directions, 
respectively 

U 

=  column  matrix  whose  elements  are  W  and  V 

w*,  w 

=  normal  deflection  and  dimensionless  normal 
deflection,  respectively 

w 

=  function  defined  by  eq.  (2-19) 

V 

=  function  defined  by  eq.  (2-20) 

V 

=  Poisson's  ratio  (take  v  =  0.3) 

T 

=  dimensionless  time 

P 

=  mass  density  of  the  shell 

X 

=  aspect  ratio  a/b 

AC,  An 

=  spacing  dimensions 

At 

=  time  increment 

a*,  a*,  t* 
x*  y’  xy 

=  membrane  stresses  in  the  middle  surface 

°1  *  °2»  t12 

=  dimensionless  membrane  in  the  middle  surface 

ar*  ae 

=  membrane  stresses  in  the  middle  surface 

as’  °e 

=  dimensionless  membrane  stresses 

Bk 

=  independent  random  number 

C1 »  c2’  c3 

=  constants 

D  (z) 

=  parabolic  cylinder  function 

e 

=  difference  between  a  nonlinear  system  and  its 
equivalent  linear  system 

E[  1 

=  expectation  operator 

Go 

f(t) 

=  power  spectral  density  of  the  load 

=  generalized  coordinate 

Ah 

=  basic  time  increment 

p(y-[»  y2) 

=  a  joint  probability  density  for  y-j  and 

yr  y2 

-  dynamical  variables 

vk 

=  a  white  Gaussian  number 

a>  B,  Y »  w0 

=  constants 

C 

=  fraction  of  critical  damping 

(0 

=  frequency 

2 

U) 

oe 

=  equivalent  linear  stiffness 

At 

=  integration  time  increment 

2 

°w 

2 

°f 

2 

°fo 

'f'(C.n) ,  *(s) 

=  variance  of  the  deflection 

=  variance  of  the  generalized  coordinate  f 

=  variance  of  f  in  a  linear  system 

=  coordinate  functions 
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Figure  2-Zt  Geometry  of  Shallow  Spherical  Shell 
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Figure  2~3:  Network  for  Rectangular  Plate 


Figure  2-4:  Descrete  Points  at  Time 
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Figure  2-6 »  Deflection  of  Clamped  Spherical  Shell 
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Figure  ~-l»  Stationary  '-ean-Square  Response  of  Simply  Supported  Square  Flats  to 
White  Noise 
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Dimensionless  Power  Spectral  Density  G0 

Figure  4-2:  Stationary  Mean-Square  Response  of  Simply  Supported  Rectangular  Plate 
(aspect  ratio  2)  to  White  Noise 
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Figure  4-3 »  Stationary  Mean-square  Response  of  Clamped  Square  Plate  to  White  Noise 
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Figure  4-6  s  Stationary  Wean-square  Response  of  Clamped  Circular  Plate  to 


Simulation  by  the  Runge-Kutta 
Method 
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Dimensionless  Power  Spectral  Density  G0 
Figure  4-7*  Stationary  Variance  of  Simply  Supported  Cylindrical  Shell  (Aspect 
ratio  /\=1»  k,  =  0  and  k«=5)  to  White  Noise 
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Fokker-Planck  Approach 
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The  Compatibility  Equation  at  the  Center  of  the  Shell  s=0 


The  compatibility  equation  in  dimensionless  form  is 

A  +  2A  .  +  -je. — A  -m  .k(B_  +  A 

3s4  sasJ  s  as^  s  as  as^  sas  sas  as^ 


(A—  1 ) 


At  s=0. 


ii i  -  0 

as 


(A-2) 


Therefore,  using  L'Hospital's  rule, 
nim  aw  _  a2w 

slV  nr  - 

The  right  hand  side  of  eq.(A-l)  becomes 

R.H.S.  =  -  (-^-y)  -2k(^-y) 

as^  as^ 

which  is  bounded  at  s=0.  Expanding  F  into  a  Taylor  series  at  s=0,  we  have 


(A-3) 


(A-4) 


a0  +  ais  +  a2s 


’+  ans 


(A-  5 ) 


where 


a0  ■  Fs=0 
a,  ■  (3F/3s)s=0 

.  (A-6; 

an  =  n!^n,:/asn)s=0 

Substituting  eq.(A-5)  into  the  left  hand  side  of  eq.  (A-l),  v^F  is  expressed 
v^F  =  a-|S-3  +  9a3s~^  +  64a^  +  165agS  +  •••  (A-  7 ) 

Since  the  right  hand  side  of  eq.(A-l)  is  bounded  at  s=0,  the  left  hand  side 
must  be  bounded.  From  eq.(A-7)  we  must  have 

a]  =  a3  =  0  (A-8) 

Therefore, 

1  imv4F  =  64a,  = 

s  0  asH  s=0  (A-9) 

The  compatibility  equation  at  s=0  is  expressed  as 

I  A-  -  -  (A)  -  2k(-^|) 

asH  as  as^ 


( A- 1 0 ) 
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